# libraries
library(readr)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(gganimate)
library(gifski)
library(dplyr)
library(tidyr)
library(ggthemes)
library(ggpubr)
library(lme4)
library(ggeffects)
library(lmerTest)
library(car)
library(ggplot2)
library(survival)
library(coxme)
library(survminer)
library(ggfortify)
library(broom)
library(ehahelper)
library(dena)
library(tidyr)
sessionInfo()
R version 4.3.0 (2023-04-21)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/London
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] dena_0.1.0 ehahelper_0.3.9999 broom_1.0.4 ggfortify_0.4.16
[5] survminer_0.4.9 coxme_2.2-18.1 bdsmatrix_1.3-6 survival_3.5-5
[9] car_3.1-2 carData_3.0-5 lmerTest_3.1-3 ggeffects_1.2.2
[13] lme4_1.1-33 Matrix_1.5-4 ggpubr_0.6.0 ggthemes_4.2.4
[17] tidyr_1.3.0 gifski_1.6.6-1 gganimate_1.0.8 ggplot2_3.4.2
[21] sjlabelled_1.2.0 sjmisc_2.8.9 sjPlot_2.8.14 readr_2.1.4
[25] data.table_1.14.8 dplyr_1.1.2
loaded via a namespace (and not attached):
[1] gridExtra_2.3 rlang_1.1.1 magrittr_2.0.3 compiler_4.3.0
[5] mgcv_1.8-42 vctrs_0.6.2 stringr_1.5.0 pkgconfig_2.0.3
[9] crayon_1.5.2 fastmap_1.1.1 backports_1.4.1 labeling_0.4.2
[13] KMsurv_0.1-5 effectsize_0.8.3 utf8_1.2.3 rmarkdown_2.21
[17] tzdb_0.3.0 haven_2.5.2 nloptr_2.0.3 purrr_1.0.1
[21] bit_4.0.5 xfun_0.39 cachem_1.0.8 jsonlite_1.8.4
[25] progress_1.2.2 tweenr_2.0.2 parallel_4.3.0 prettyunits_1.1.1
[29] R6_2.5.1 bslib_0.4.2 stringi_1.7.12 pkgload_1.3.2
[33] boot_1.3-28.1 jquerylib_0.1.4 numDeriv_2016.8-1.1 estimability_1.4.1
[37] Rcpp_1.0.10 knitr_1.42 modelr_0.1.11 zoo_1.8-12
[41] parameters_0.21.0 splines_4.3.0 tidyselect_1.2.0 yaml_2.3.7
[45] rstudioapi_0.14 abind_1.4-5 lattice_0.21-8 tibble_3.2.1
[49] withr_2.5.0 bayestestR_0.13.1 evaluate_0.20 survMisc_0.5.6
[53] pillar_1.9.0 insight_0.19.1 plotly_4.10.1 generics_0.1.3
[57] vroom_1.6.3 hms_1.1.3 munsell_0.5.0 scales_1.2.1
[61] minqa_1.2.5 xtable_1.8-4 glue_1.6.2 emmeans_1.8.5
[65] lazyeval_0.2.2 tools_4.3.0 ggsignif_0.6.4 forcats_1.0.0
[69] mvtnorm_1.1-3 cowplot_1.1.1 grid_4.3.0 datawizard_0.7.1
[73] colorspace_2.1-0 nlme_3.1-162 eha_2.10.3 performance_0.10.3
[77] cli_3.6.1 km.ci_0.5-6 fansi_1.0.4 viridisLite_0.4.2
[81] sjstats_0.18.2 gtable_0.3.3 rstatix_0.7.2 sass_0.4.6
[85] digest_0.6.31 htmlwidgets_1.6.2 farver_2.1.1 htmltools_0.5.5
[89] lifecycle_1.0.3 httr_1.4.5 bit64_4.0.5 MASS_7.3-60
# read the data file "ADD_I"
ADD_I <- read_csv("Tables/ADD_I.csv", col_types = cols(Island = col_factor(levels = c("CN","AR", "CE", "DS", "FR"))))
# create the earliest born individuals on each island
Isfy <- ADD_I %>% select(founderyear,Island) %>% filter(!duplicated(Island))
# create a data frame with means and sd of each morphometric
ADD_mean <- ADD_I %>% group_by(Island,BirthYear) %>% summarise(RTsd = sd(RightTarsus,na.rm=TRUE),WLsd = sd(WingLength,na.rm=TRUE),BMsd = sd(BodyMass,na.rm=TRUE),HBsd=sd(HeadBill,na.rm=TRUE),RightTarsus=mean(RightTarsus,na.rm=TRUE),WingLength=mean(WingLength,na.rm=TRUE),BodyMass=mean(BodyMass,na.rm=TRUE),HeadBill=mean(HeadBill,na.rm=TRUE))
`summarise()` has grouped output by 'Island'. You can override using the `.groups` argument.
####RightTarsus ----
IRTlmer <- lmer(RightTarsus ~ Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IRTlmer)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: RightTarsus ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +
(1 | BirdID)
Data: ADD_I
REML criterion at convergence: 10511
Scaled residuals:
Min 1Q Median 3Q Max
-5.0560 -0.4445 0.0052 0.4324 5.5726
Random effects:
Groups Name Variance Std.Dev.
BirdID (Intercept) 0.3048654 0.55215
Observer (Intercept) 0.0562847 0.23724
Ageclass (Intercept) 0.0009875 0.03142
Residual 0.1407944 0.37523
Number of obs: 6154, groups: BirdID, 3258; Observer, 82; Ageclass, 5
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.417e+01 6.083e-02 2.713e+02 397.385 < 2e-16 ***
IslandAR -7.215e-02 1.686e-01 3.839e+03 -0.428 0.668783
IslandCE -2.510e-01 9.829e-02 3.925e+03 -2.553 0.010704 *
IslandDS 2.519e-01 8.080e-02 3.408e+03 3.117 0.001842 **
IslandFR 1.869e-01 1.236e-01 2.500e+03 1.512 0.130545
fYear 8.838e-03 2.530e-03 2.174e+03 3.493 0.000488 ***
SexMale 1.596e+00 2.244e-02 3.211e+03 71.130 < 2e-16 ***
IslandAR:fYear -6.760e-04 1.024e-02 3.539e+03 -0.066 0.947388
IslandCE:fYear 1.764e-03 6.196e-03 4.079e+03 0.285 0.775832
IslandDS:fYear -2.910e-03 5.840e-03 4.085e+03 -0.498 0.618312
IslandFR:fYear 1.554e-03 1.756e-02 1.981e+03 0.089 0.929486
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) IslnAR IslnCE IslnDS IslnFR fYear SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR -0.163
IslandCE -0.309 0.111
IslandDS -0.473 0.113 0.218
IslandFR -0.331 0.081 0.150 0.225
fYear -0.779 0.188 0.339 0.504 0.366
SexMale -0.177 -0.029 0.008 -0.033 -0.024 -0.016
IslndAR:fYr 0.149 -0.922 -0.096 -0.101 -0.076 -0.200 0.026
IslndCE:fYr 0.261 -0.094 -0.915 -0.173 -0.126 -0.323 -0.010 0.097
IslndDS:fYr 0.314 -0.065 -0.138 -0.849 -0.208 -0.379 0.008 0.069 0.123
IslndFR:fYr 0.102 -0.025 -0.050 -0.108 -0.850 -0.143 0.020 0.029 0.048 0.215
#predicting the model
IRTlmerdata <- ggpredict(IRTlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,RightTarsus=predicted,Island=group,Sex=facet)
IRTlmerdata2 <- merge(IRTlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear) %>% summarise(RightTarsus = mean(RightTarsus, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))
#plotting the model (shown as ggarrange below)
IRTmod <- ggplot(IRTlmerdata2, aes(x = BirthYear, y = RightTarsus, fill = Island)) +
geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=RightTarsus,colour=Island), size = 2) +
geom_errorbar(data=ADD_mean,aes(ymin=RightTarsus-RTsd,ymax=RightTarsus+RTsd, colour=Island),alpha=0.7) +
geom_smooth(aes(colour = Island),method = "loess", se = FALSE) +
geom_ribbon(aes(ymin=RightTarsus-std.error,ymax=RightTarsus+std.error),alpha=0.5) +
labs(x = "Birth Year", y= "Tarsus Length (mm)") +
xlim(1990,2022) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
#### Body Mass ----
IBMlmer <- lmer(BodyMass ~ Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IBMlmer)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: BodyMass ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +
(1 | BirdID)
Data: ADD_I
REML criterion at convergence: 16677.7
Scaled residuals:
Min 1Q Median 3Q Max
-3.6968 -0.5705 -0.0271 0.5440 4.2625
Random effects:
Groups Name Variance Std.Dev.
BirdID (Intercept) 0.24508 0.4951
Observer (Intercept) 0.05818 0.2412
Ageclass (Intercept) 0.14657 0.3828
Residual 0.53384 0.7306
Number of obs: 6554, groups: BirdID, 3398; Observer, 91; Ageclass, 5
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.438e+01 1.846e-01 5.339e+00 77.909 2.29e-09 ***
IslandAR 1.189e+00 1.570e-01 3.120e+03 7.577 4.63e-14 ***
IslandCE 2.792e-01 1.214e-01 3.666e+03 2.299 0.02158 *
IslandDS 7.385e-01 1.037e-01 1.096e+03 7.122 1.92e-12 ***
IslandFR 5.264e-01 1.668e-01 1.471e+03 3.156 0.00163 **
fYear -2.623e-03 3.107e-03 5.626e+02 -0.844 0.39893
SexMale 1.433e+00 2.635e-02 2.907e+03 54.405 < 2e-16 ***
IslandAR:fYear -7.870e-02 1.043e-02 2.346e+03 -7.546 6.38e-14 ***
IslandCE:fYear -3.708e-02 7.935e-03 4.479e+03 -4.673 3.05e-06 ***
IslandDS:fYear -5.069e-02 7.678e-03 2.097e+03 -6.601 5.14e-11 ***
IslandFR:fYear -5.400e-02 2.379e-02 1.097e+03 -2.270 0.02343 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) IslnAR IslnCE IslnDS IslnFR fYear SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR -0.071
IslandCE -0.107 0.098
IslandDS -0.176 0.107 0.184
IslandFR -0.113 0.071 0.120 0.187
fYear -0.307 0.188 0.294 0.462 0.320
SexMale -0.068 -0.016 0.007 -0.035 -0.027 -0.017
IslndAR:fYr 0.065 -0.869 -0.082 -0.095 -0.067 -0.205 0.016
IslndCE:fYr 0.092 -0.086 -0.910 -0.144 -0.097 -0.283 -0.007 0.092
IslndDS:fYr 0.117 -0.064 -0.114 -0.847 -0.181 -0.345 0.015 0.069 0.098
IslndFR:fYr 0.035 -0.024 -0.039 -0.097 -0.868 -0.127 0.023 0.027 0.035 0.201
#predicting the model
IBMlmerdata <- ggpredict(IBMlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,BodyMass=predicted,Island=group,Sex=facet)
IBMlmerdata2 <- merge(IBMlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear) %>% summarise(BodyMass = mean(BodyMass, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))
#plotting the model
IBMmod <- ggplot(IBMlmerdata2, aes(x = BirthYear, y = BodyMass, fill = Island)) +
geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=BodyMass,colour=Island), size = 2) +
geom_errorbar(data=ADD_mean,aes(ymin=BodyMass-BMsd,ymax=BodyMass+BMsd, colour=Island),alpha=0.7) +
geom_smooth(aes(colour = Island),method = "lm", se = FALSE) +
geom_ribbon(aes(ymin=BodyMass-std.error,ymax=BodyMass+std.error),alpha=0.5) +
labs(x = "Birth Year", y= "Body Mass (g)") +
xlim(1990,2022) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
#### Wing Length ----
IWLlmer <- lmer(WingLength ~ Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
Warning: Model failed to converge with max|grad| = 0.0760292 (tol = 0.002, component 1)
summary(IWLlmer)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: WingLength ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +
(1 | BirdID)
Data: ADD_I
REML criterion at convergence: 23050.6
Scaled residuals:
Min 1Q Median 3Q Max
-3.8213 -0.5350 0.0161 0.5335 4.0139
Random effects:
Groups Name Variance Std.Dev.
BirdID (Intercept) 0.95940 0.9795
Observer (Intercept) 0.40157 0.6337
Ageclass (Intercept) 0.08901 0.2983
Residual 1.26489 1.1247
Number of obs: 6510, groups: BirdID, 3408; Observer, 91; Ageclass, 5
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 6.592e+01 1.905e-01 1.616e+01 346.012 < 2e-16 ***
IslandAR 2.490e-01 2.795e-01 3.717e+03 0.891 0.37296
IslandCE -1.869e-01 2.098e-01 4.116e+03 -0.891 0.37310
IslandDS -5.162e-01 1.865e-01 2.916e+03 -2.767 0.00569 **
IslandFR 1.871e-01 2.876e-01 2.198e+03 0.650 0.51549
fYear -2.900e-02 5.796e-03 1.679e+03 -5.004 6.22e-07 ***
SexMale 2.589e+00 4.624e-02 3.176e+03 55.991 < 2e-16 ***
IslandAR:fYear -4.495e-02 1.871e-02 3.681e+03 -2.403 0.01632 *
IslandCE:fYear 2.268e-03 1.354e-02 4.598e+03 0.168 0.86696
IslandDS:fYear 4.375e-02 1.350e-02 4.156e+03 3.241 0.00120 **
IslandFR:fYear 2.893e-02 4.134e-02 2.172e+03 0.700 0.48407
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) IslnAR IslnCE IslnDS IslnFR fYear SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR -0.123
IslandCE -0.198 0.099
IslandDS -0.315 0.101 0.190
IslandFR -0.225 0.072 0.131 0.197
fYear -0.551 0.177 0.302 0.464 0.344
SexMale -0.118 -0.015 0.005 -0.033 -0.026 -0.015
IslndAR:fYr 0.115 -0.873 -0.083 -0.092 -0.070 -0.200 0.013
IslndCE:fYr 0.174 -0.088 -0.910 -0.152 -0.112 -0.299 -0.007 0.096
IslndDS:fYr 0.218 -0.061 -0.121 -0.851 -0.195 -0.359 0.012 0.068 0.107
IslndFR:fYr 0.074 -0.023 -0.043 -0.099 -0.857 -0.135 0.024 0.028 0.040 0.217
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0760292 (tol = 0.002, component 1)
#predicting the model
IWLlmerdata <- ggpredict(IWLlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,WingLength=predicted,Island=group,Sex=facet)
IWLlmerdata2 <- merge(IWLlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11)) %>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear)%>% summarise(WingLength = mean(WingLength, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))
#plotting the model
IWLmod <- ggplot(IWLlmerdata2, aes(x = BirthYear, y = WingLength, fill = Island)) +
geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=WingLength,colour=Island), size = 2) +
geom_errorbar(data=ADD_mean,aes(ymin=WingLength-WLsd,ymax=WingLength+WLsd, colour=Island),alpha=0.7) +
geom_smooth(aes(colour = Island),method = "lm", se = FALSE) +
geom_ribbon(aes(ymin=WingLength-std.error,ymax=WingLength+std.error),alpha=0.5) +
labs(x = "Birth Year", y= "Wing Length (mm)") +
xlim(1990,2022) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
#### Head + Bill ----
IHBlmer <- lmer(HeadBill ~ Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IHBlmer)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: HeadBill ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +
(1 | BirdID)
Data: ADD_I
REML criterion at convergence: 9601.4
Scaled residuals:
Min 1Q Median 3Q Max
-4.1841 -0.4234 0.0279 0.4642 5.4733
Random effects:
Groups Name Variance Std.Dev.
BirdID (Intercept) 0.22985 0.4794
Observer (Intercept) 0.04913 0.2217
Ageclass (Intercept) 0.04152 0.2038
Residual 0.14507 0.3809
Number of obs: 5937, groups: BirdID, 3120; Observer, 80; Ageclass, 5
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.588e+01 1.085e-01 7.815e+00 330.737 <2e-16 ***
IslandAR 4.153e-01 2.063e-01 1.204e+03 2.013 0.0443 *
IslandCE -1.041e-01 9.292e-02 3.741e+03 -1.121 0.2624
IslandDS 8.141e-02 7.652e-02 3.252e+03 1.064 0.2874
IslandFR -7.000e-02 1.160e-01 2.381e+03 -0.603 0.5464
fYear -1.237e-03 2.492e-03 2.573e+03 -0.496 0.6196
SexMale 8.420e-01 2.075e-02 3.016e+03 40.579 <2e-16 ***
IslandAR:fYear -2.977e-02 1.559e-02 5.975e+02 -1.910 0.0566 .
IslandCE:fYear -8.098e-03 5.841e-03 3.978e+03 -1.386 0.1657
IslandDS:fYear 1.404e-03 5.456e-03 3.918e+03 0.257 0.7969
IslandFR:fYear 1.196e-02 1.638e-02 1.925e+03 0.730 0.4655
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) IslnAR IslnCE IslnDS IslnFR fYear SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR -0.066
IslandCE -0.177 0.079
IslandDS -0.277 0.088 0.233
IslandFR -0.197 0.067 0.163 0.243
fYear -0.444 0.155 0.351 0.525 0.386
SexMale -0.089 -0.017 -0.001 -0.036 -0.029 -0.021
IslndAR:fYr 0.046 -0.944 -0.049 -0.062 -0.051 -0.133 0.009
IslndCE:fYr 0.152 -0.064 -0.918 -0.188 -0.138 -0.333 -0.002 0.047
IslndDS:fYr 0.187 -0.049 -0.147 -0.851 -0.220 -0.393 0.012 0.041 0.131
IslndFR:fYr 0.062 -0.021 -0.054 -0.113 -0.848 -0.148 0.023 0.019 0.050 0.218
#predicting the model
IHBlmerdata <- ggpredict(IHBlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,HeadBill=predicted,Island=group,SexEstimate=facet)
IHBlmerdata2 <- merge(IHBlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 22)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear)%>% group_by(Island,BirthYear)%>% summarise(HeadBill = mean(HeadBill, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))
#plotting the model
IHBmod <- ggplot(IHBlmerdata2, aes(x = BirthYear, y = HeadBill, fill = Island)) +
geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=HeadBill,colour=Island), size = 2) +
geom_errorbar(data=ADD_mean,aes(ymin=HeadBill-HBsd,ymax=HeadBill+HBsd, colour=Island),alpha=0.7) +
geom_smooth(aes(colour = Island),method = "lm", se = FALSE) +
geom_ribbon(aes(ymin=HeadBill-std.error,ymax=HeadBill+std.error),alpha=0.25) +
labs(x = "Birth Year", y= "Head + Bill Length (mm)") +
xlim(1990,2022) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
Imods <- ggarrange(IRTmod,IBMmod,IWLmod,IHBmod, common.legend=TRUE, labels = c("A","B","C","D"))
`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'
Imods
ADD_R <- read_csv("Tables/ADD_R.csv", col_types = cols(SexEstimate = col_factor()))
ADD_Real <- as.data.frame(ADD_R)
ADD_Real <- ADD_Real %>% select(RightTarsus,BodyMass,WingLength,HeadBill,BirdID,SexEstimate,BirthYear,Lifespan,SurvStatus,ReproductiveOutput)
ADD_tidy <- ADD_Real %>% filter(!is.na(RightTarsus) & !is.na(BodyMass) & !is.na(WingLength) & !is.na(HeadBill) & !is.na(SexEstimate) & !is.na(Lifespan) ) %>% mutate(BodyCondition = BodyMass/RightTarsus) %>% mutate(SexEstimate = as.factor(SexEstimate))
###Lifespan analysis with coxme----
lscox0 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ BodyMass + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox1 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ poly(BodyMass,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox2 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus*poly(BodyMass,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox3 <- coxme(Surv(Lifespan,SurvStatus) ~ poly(BodyCondition,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox4 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ poly(BodyMass,2)*SexEstimate + WingLength + HeadBill + (1|BirthYear), data= ADD_tidy)
AIC(lscox0,lscox1,lscox2,lscox3,lscox4) #quadratic body mass is best
BIC(lscox0,lscox1,lscox2,lscox3,lscox4)
# lscox4 was selected based on the best AIC and BIC values
####predict lifespan----
ADD_round <- ADD_R %>%mutate(across(c(RightTarsus, BodyMass, HeadBill, WingLength), round, 0)) %>% distinct(RightTarsus, BodyMass, HeadBill, WingLength,SexEstimate) %>% filter(!is.na(RightTarsus) & !is.na(BodyMass) & !is.na(WingLength) & !is.na(HeadBill) & !is.na(SexEstimate))
ADD_end <- bind_rows(ADD_tidy,ADD_round)
rs <- predict_coxme(lscox4,newdata = ADD_end,type="risk", se.fit=TRUE)
rsend <- ADD_end %>% add_columns(rr = rs$fit[,1], se.fit = rs$se.fit[,1])
ADD_tail <- rsend %>% tail(nrow(ADD_round)) %>% mutate(SexEstimate = as.factor(SexEstimate)) %>% mutate(low = exp(log(rr) - 1.96 * se.fit),high = exp(log(rr) + 1.96 * se.fit)) %>% group_by(BodyMass,SexEstimate) %>% summarise(mean_rr=mean(rr),mean_se.fit=mean(se.fit), mean_low=mean(low),mean_high=mean(high))
`summarise()` has grouped output by 'BodyMass'. You can override using the `.groups` argument.
coxmeplot <- ggplot(ADD_tail,aes(x=BodyMass,y=mean_rr)) +
geom_line(aes(colour=SexEstimate)) +
geom_ribbon(aes(fill=SexEstimate,ymin=mean_rr-mean_se.fit,ymax=mean_rr+mean_se.fit), alpha = 0.25) +
labs(x = "Body Mass (g)", y = "Mortality Risk Score") +
guides(fill=guide_legend(title="Sex"), colour=guide_legend(title="Sex")) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind(labels = c("Female", "Male")) +
scale_fill_colorblind(labels = c("Female", "Male"))
ROglmF <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill ,data = ADD_tidy, family = "poisson")
ROglmR <- glm(ReproductiveOutput ~ poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm0 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill + Lifespan + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm1 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm2 <- glm(ReproductiveOutput ~ poly(RightTarsus,2)+ BodyMass + WingLength + HeadBill + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm3 <- glm(ReproductiveOutput ~ RightTarsus+ poly(BodyMass,2) + WingLength + HeadBill + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm4 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + poly(WingLength,2) + HeadBill + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm5 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + poly(HeadBill,2) + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
AIC(ROglmF,ROglmR,ROglm0,ROglm1,ROglm2,ROglm3,ROglm4,ROglm5) # 4 best, 1 second
BIC(ROglmF,ROglmR,ROglm0,ROglm1,ROglm2,ROglm3,ROglm4,ROglm5) # 1 best
vif(ROglm1)
GVIF Df GVIF^(1/(2*Df))
RightTarsus 3.252977 1 1.803601
BodyMass 2.553710 1 1.598033
WingLength 2.550319 1 1.596972
HeadBill 2.274822 1 1.508251
poly(Lifespan, 2) 1.212419 2 1.049333
SexEstimate 3.638709 1 1.907540
BirthYear 1.266770 1 1.125509
summary(ROglm1)
Call:
glm(formula = ReproductiveOutput ~ RightTarsus + BodyMass + WingLength +
HeadBill + poly(Lifespan, 2) + SexEstimate + BirthYear, family = "poisson",
data = ADD_tidy)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 55.943199 6.971746 8.024 1.02e-15 ***
RightTarsus -0.004860 0.032759 -0.148 0.88207
BodyMass 0.075016 0.028836 2.602 0.00928 **
WingLength 0.015426 0.015372 1.004 0.31561
HeadBill -0.002388 0.038924 -0.061 0.95108
poly(Lifespan, 2)1 23.108050 1.113788 20.747 < 2e-16 ***
poly(Lifespan, 2)2 -5.454150 0.702063 -7.769 7.93e-15 ***
SexEstimate1 -0.042798 0.066926 -0.639 0.52251
BirthYear -0.028461 0.003275 -8.692 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 2018.43 on 736 degrees of freedom
Residual deviance: 879.58 on 728 degrees of freedom
(629 observations deleted due to missingness)
AIC: 3146.9
Number of Fisher Scoring iterations: 5
ROplotdata <- ADD_tidy %>% mutate(BodyMass = case_when(BodyMass > 16.5 ~ "16.5", BodyMass < 14.5 ~ "14.5",BodyMass < 16.5 & BodyMass > 14.5 ~ "15.5")) %>% filter(!is.na(BodyMass))
ROglmdat <- ggpredict(ROglm1, terms = c("Lifespan","BodyMass")) %>% rename(Lifespan=x,ReproductiveOutput=predicted,BodyMass=group)
ROglmplot <- ggplot(ROglmdat, aes(x=Lifespan,y=ReproductiveOutput)) +
geom_point(data = ROplotdata, mapping=aes(x=Lifespan,y=ReproductiveOutput,colour=BodyMass),position="jitter") +
geom_smooth(aes(colour=BodyMass), se = FALSE) +
geom_ribbon(aes(ymin=conf.low,ymax=conf.high,fill=BodyMass), alpha = 0.25) +
labs(x = "Lifespan (years)", y= "Total Number of Offspring") +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
fitplot <- ggarrange(coxmeplot,ROglmplot, labels = c("A","B"))
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
fitplot
Thank you for going through the code and supporting open research.